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Abstract. The evolution of dynamical perturbations is exam- 
ined in nuclear multifragmentation in the frame of Vlasov equation. 
Both plane wave and bubble type of perturbations are investigated 
in the presence of surface (Yukawa) forces. An energy condition is 
given for the allowed type of instabilities and the time scale of the 
exponential growth of the instabilities is calculated. The results are 
compared to the mechanical spinodal region predictions. 
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1. Introduction 



Let us consider a spherically expanding nuclear system in the metastable nuclear 
fluid phase when it reaches the freeze-out at time r/y. Although at the freeze- 
out the fermionic degrees of freedom are frozen-out, and internucleon collisions 
cease, softer long range nuclear interactions are still effective, and represented by 
a nuclear mean field potential, U(f). 

We will assume that the system, both before and after freeze-out undergoes 
a spherical, scaling expansion. Such an expansion can be represented by a four- 
velocity field, = x^/t, where r = y/t 2 — x 2 — y 2 — z 2 , for the internal regions 
of our expanding system (but not for the external surface). This flow pattern is 
invariant under Lorentz transformation, i.e., the points of the interior of our ex- 
panding system are physically identical and indistinguishable from one another. 
Consequently in the interior, in the Local Rest (LR) frame all thermodynamical 
and fluid-dynamical quantities are equal, and from the point of view of instabil- 
ities all internal points are equivalent. We will exploit these symmetry features 
although in the calculations we will use a non-relativistic approximation. 

Let us assume that the nucleon phase space distribution before, and at the 
freeze-out, Tf r , is a Fermi distribution: 

fo(Tf r ,r,p) = cjl + exp 

where this form assumes a correlation between the momentum and radial dis- 
tribution arising from radial expansion, C = g/(2irh) 3 is the normalization, and 
g is the degeneracy of nucleons, so that the proper (LR) density is n (rj r .,r) = 
/ d 3 p fo{Tf r , r,p). We assume that in the interior of the collision zone the freeze- 
out density is constant: n (rf r , r) — nyy. In the center of the collision zone this is 
an ideal Fermi distribution, while at finite radii, |r*|, the distribution is boosted 
(using non-relativistic, Galilei transformation), with a radially directed and ra- 
dially linearly increasing flow velocity of v — f/r. We can also introduce the LR 
momentum: P(r, r ) = p — mf/r. 

Furthermore, let us assume that after the freeze-out for r > Tj r , the system 
expands homogeneously according to the collisionless Vlasov equation. The dis- 
tribution function, /, is the solution of the Vlasov equation with a mean-field 
potential, U(r), 

dr m df df dp 
In the special case of a homogeneous system, where the last term vanishes, for 
such a free coasting expansion in the local rest frame is just obtained by replacing 
p by rp/r fr in Eq. © § §: 

fo(r,r,p) = c\l+exp 
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where T e // = T/ r (r/ r /r) 2 and // e // = /i/ r (r/, r /r) 2 . The condition, that the ratio of 
the chemical potential and the temperature is constant during the expansion, in a 
usual thermodynamical system, corresponds to an adiabatic process. The density 
of the system changes with time in this inertial expansion as uq(t) = n/ r (r/ r /r) 3 . 
This solution is valid starting from the freeze-out, T/ r , and until inhomogeneities 
will spontaneously develop in the system at some threshold time, r^. Before 
this time small perturbations will smooth out due to the mean field potential, 
while after this threshold time the density dependent mean field will enhance 
fluctuations. So after this threshold time the density will not be homogeneous 
any more. 

Note that this post freeze-out distribution, /, is not a thermal equilibrium 
distribution function, and the effective parameters, T e ff and /i e // are just carrying 
the memory of the last equilibrium thermal parameters, Tf r and fif r , but these 
are not the usual thermodynamical parameters. This can be easily seen if the 
expansion is not spherically symmetric [|TJ. 

If we would have a thermal expansion following rj r , the Equation of State 
(EOS) would determine the time dependence of the physical temperature in an 
adiabatic expansion. Generally this would not coincide with T e ff = Tf r (rf r /r) 2 , 
only in the case of a large system, where the flow dominates the energy. For 
example if the EOS is that of an ideal Stefan-Boltzmann gas, (dp/de = cjj, 
Cq = 1/3 and e = cT 4 ), then T(r) = Tf r (r / Tf r )~ 3c Q , which differs from T e ff. 

We study the stability of the system and the occurrence of instabilities arising 
from the mean field. Such an instability may lead to a rapid multifragmentation 
of our system. 



2. Instabilities 



Let us consider a small perturbation in the expanding system. The amplitude of 
this perturbation may grow, decrease or oscillate depending on the conditions. 
In the presence of such a perturbation the phase space distribution is: 

f(r,r,p) = f + fi(T,f,p ), 

with the normalization n(r, f) — j d 3 p (/ + f\) and ni(r, f ' ) = / d 3 p fx, where 
the unperturbed density n is homogeneous, Vn = 0, and f\ should be a local 
spherical perturbation which is a solution of the Vlasov equation 
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The solution of the Vlasov equation is treated in details in Ref . . Here we sepa- 
rate the Vlasov equation into two equations, one for f and one for fi. Separating 
non-vanishing dominant zeroth order terms we get the equation 



dfo V dfp 
dr m dr 
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which is satisfied by fo as given in Eq. (|3|). The first order terms yield the 
linearized equation 

dr m df dp 
We intend to find perturbations which grow, leading to instabilities of the system. 
Some modes of growing perturbations may arise in thermal surrounding, and 
their rate is determined by thermal and viscous damping |3|, |6|, 0, [8[]. These 
are usually slower, and so other faster processes may come into play also. Here 
we intend to study non-thermal, growing perturbations, which may occur after 
the thermal freeze-out only, but they can be faster than the thermally damped 
processes |], |, |5j, §, [TO], [TT| . The growth rate of such perturbations is determined 
by the long range nuclear mean field potential, U(r). Usually different non- 
thermal channels of instability open only when a given time is passed after the 
freeze-out at Tf r , and we reach a threshold time, r th . 

There are two characteristic time scales in the system: (i) the longer post 
freeze-out expansion between T/ r and r th , and (ii) the rapid growth of instability 
which develops after r^. When studying the dynamics of rapidly growing insta- 
bilities we can usually neglect the much slower dynamics of the post freeze-out 
expansion. 

Different configurations can and should be taken into account when studying 
instabilities in a quenched (supercooled) system. 



2.1. Plane wave perturbation 

The stability of the Vlasov equation against plane wave perturbations was exam- 
ined in detail by different groups recently |L2|, |13| . If we expand the perturbation, 

A, as 

h(r,p,t)=J2 fk(p,t)^=^, (5) 
k vS2 

and search for the solution of fk{p,t) in the form 

fk( P ,t) = f k ,M e tuJt , 



according to Ref. [12] we get the dispersion relation uj{k) 

dU(k) r d 3 p {kp) 2 df 



dn J (2irh) 3 (kp)-m 2 to 2 de 
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where / is the static uniform solution, and U(k) denotes the Fourier component 
of the effective field U(f ). The mode corresponding to wavenumber k will be 
unstable, when the u(k) frequency becomes imaginary. It was also shown in 
Ref. , that for zero temperature the condition of instability can be written as 

2 dU(k) 
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with €f Fermi kinetic energy. The expression goes over for k — > into the 
condition of mechanical instability (i.e., the compressibility becomes negative) 

2 dU(n) 

3 on 

where the potential U(n) depends only on the homogeneous density. 

In Ref. |l3j the instability condition was examined for a one dimensionally 
expanding ground state system. The acting force was a Skyrme force with Yukawa 
surface term. The resulting instability condition reads as Eq. (|7]), where 

dU ( k ) . n« , t i 1 \ / i n\ (J 4 ^0 



-2 / 9 + T ( ( 7+l)(a + 2K- 



dn ' y n A 70 ^ + #2 ' 

and /3, 7, cr, Vo and \x are the Skyrme and Yukawa force parameters (see Eq. (|TTD ) . 

In case of a uniformly expanding system the same condition of instability for 
a plane wave perturbation can also be obtained in a fashion similar to the spher- 
ical case discussed later, see Appendix A. This approach leads to the following 
condition 



2wT " mC ^ *^*> 7, _, , (8) 

; dn J (y + s)[l + e»-'W 2 >'- 1 



where s = tuk r /(2k Tf r rj r ) , k = iuj and dU(k)/dn is the same expression as 



in Ref. 13 



2.2. Spherical drop/bubble perturbations 



In the following we want to study more realistic perturbations instead of plane 
waves. We consider local spherical bubbles, since we believe these are the first 



instabilities which start to grow [14]. In general spherical perturbations minimize 
the surface and surface energy, so these can be formed the earliest. (On the 
other hand plane wave perturbations may grow more rapidly at later stages with 
stronger driving forces due to the increased surface.) 

Consider a spherical drop with a surface density profile exponentially decreas- 
ing characterized by the parameter k, a central density, n c , and radius R(r) 



fi — fs 



n c exp 



-k^f- (r - R(t))] r > R(t) 

r < R(t) 



(9) 



where R(r) — R* y- e K<yT ~~ Tth ^ is the r dependence of the radius after the threshold 
time, r t h when the droplet becomes bigger than the critical radius and it will be 
able to grow. We are interested in the initial growth rate of the radius only, so 
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that in R(t) the exponential term can be expanded into a power series for small 
r-r th , i.e., 

R(r)^R* — [1 + K ( T - r th )}. 

Tth 

Inserting this approximation into Eq. we get for the exterior part (r > R(t)) 
of the profile: 

fi = / s n c exp (kR* — ) exp ( -k—r + kR*—n{r - r th 



Since we consider small perturbations only, where the linearization of Eq. ( ^ 
holds, this solution is valid only for a short time after the instability starts to grow. 
The dispersion relation will lead to a dynamical growth factor, k, depending both 
on k and R. 



(In the plane wave expansion of the perturbation studied in Ref. [|T2| , |13fl the 
opening of the channel of the instability was indicated when u became imagi- 
nary. Thus perturbations preceding Tth did not grow. Our approach is basically 
equivalent to their one presented here, however, we wanted to emphasize that the 
instability may grow only after Tth-) 

Since / s (r, P) has a characteristic time dependence on the slow scale (i) of 
the post freeze-out expansion, we assume that its time-derivatives are negligible 
compared to other time- derivatives of f\ (r, r , p ) corresponding to rapid inequilib- 
rium processes (ii) like exp[«(r — r^)]. Furthermore, we assume that f s depends 
on the LR momentum, P only, i.e., it does not have any other dependence on r 
other than what is included in P. 

We search for such solutions of the Vlasov equation, /i(r, f,p) and f s (r,P). 
We assume that such solutions can be obtained only some time, r, after the 
freeze-out at r/ r , i.e. at Tth > T fr- Before this time thermal processes and 
thermal damping is dominant which generally lead to slower nucleation than 
post-freeze-out processes driven by the background fields. 

We can calculate the critical droplet radius, R* rit . Droplets smaller than R* Tit 
tend to disappear, while droplets larger than R* Tit may start to grow. Thus we 
will study the growth rate of critical size droplets. The critical radius, R* rit is 
calculated in Appendix E. 

The critical radius, R* rit , should be evaluated when the channel of instability 
opens at T th , and the critical droplet just starts to grow. In the 3-dimensional 
scaling expansion the critical radius scales with the overall scaling, which leads 
to a quasi-static critical radius of R* rit T/T t h- 

The critical radius, R* lit , depends on the background nucleon density at the 
opening of the instability, no (rth) = nf r (Tf r /T t h) 3 , or consequently on T t t, fur- 
thermore on the surface parameter, k, on the central density, n c , and on the 
parameters of the interaction potential. The total energy of the system can be 
simultaneously minimized by varying k and n c , and searching for an extremum 
as a function of R*. 
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As we mentioned the time-scale of the expansion is assumed to be slow com- 
pared to the time-scale of the instability, so R* can be considered as a time 
independent constant when studying the growth of instability 

The perturbation (H) satisfies the Vlasov equation both for the exterior and 
interior region, if an averaged Yukawa potential is used (see Appendix C). 

The dispersion relation for such a spherical perturbation can be written as 



(see Appendix B), where S = tuk 2 (R*) 2 t 4 / (2T fr rf r rl) . 

It is easy to see, that for S = without the surface term this condition is 
equivalent to the isothermal mechanical instability even for T ^ (Appendix 
D), that is the region of dynamical instability and of mechanical one coincide. 
Although the direct /c-dependence drops out of the dispersion relation, but since 
R* depends on k so the dispersion relation is still applicable. It is interesting 
to mention that both dispersion relations, (|S[|T0l), yield the same condition for 
evaluating the threshold time for the perturbation which is just on the boundary 
to be able to grow, i.e., s — S — 0. 

The features of the potential and the Yukawa term in it are vital in deter- 
mining the properties of the static, critical droplet or bubble. Physically we can 
consider two situations in the course of final multifragmentation. 

Depending on the beam energy, after the initial compression we reach the 
most compressed and heated up state with a definite specific entropy. This stage 
is then followed by an expansion, which is adiabatic to a good approximation. If 
the final specific entropy is smaller than the critical entropy of the nuclear liquid- 
gas phase transition, the expansion will lead to a stretched (or quenched) liquid 
state, with density no below the normal nuclear density, n^. The instabilities 
will lead then to bubble formation with an interior nuclear gas phase density, 
rii + n « 0.1 — 0An N . If on the other hand the final specific entropy exceeds the 
critical entropy, the expansion will lead to a oversaturated (or quenched) nuclear 
vapor state, with density 0.1 — 0An , below the critical nuclear density. The 
instabilities will lead now to the condensation of a nuclear liquid droplet with an 
interior nuclear density, n\ + n ~ tin- 

3. Condition for the instabilities to grow 

Equations (|8|,|iOD determine the condition for k becoming real, but its sign remains 
undefined. From equations (|5|||) it is easy to see, that the amplitude of the 
perturbation will depend on the sign of k: positive k-s will cause exponentially 
increasing perturbation, while perturbations corresponding to negative k will be 
damped rapidly. To determine the sign of the k we have to see, how the energy 




(10) 
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P (MeV fm 3 ) 


7 (MeV fm 3 ( CT+1 )) 


o 


V (MeV fm) 


/! (fill X ) 


SOFT 


1051.76 


1107.41 


1/6 


83.5 


2 


HARD 


365.0 


808.65 


1 


83.5 


2 



Table 1: The parameterization of the SOFT and HARD EOS 

changes due to the perturbation. If the configuration with the perturbation 
acquires smaller total nuclear energy (the total energy without the flow) than the 
unperturbed system can we speak about growing instabilities, that is flow can 
develop to take extra matter into the perturbation. 

In the following we consider density dependent Skyrme forces with an averaged 
Yukawa term. The total nuclear energy of the system can be written as 



E = E kin -(3 J £rn\r)+ 1 j d 3 r n CT+2 ( 



r — 



e 



Vq / d i r 1 d i r 2 n(fi)n(f 2 )— — , (11) 

ri — r 2 



where the values of (3, 7, a, Vq and \i are the same as in Ref. |HJ and summarized 
in Table 1. Let us consider a system which is initially homogeneous and has 
constant density n (r) = nj r (rj r /r) 3 . Introducing now a small perturbation in 
the density in a way that the total mass number has to be conserved, we obtain 
a density distribution 

r 

n(r,r ) = n (r) +ni(r,f ) - — (12) 

where T = J d 3 r ra 1 (r, r), Q is the volume of the system, and n\ is assumed to 
be small compared to n . If the initial configuration is such that the formation 
of a perturbation may lead to a decrease of the energy such perturbations will 
appear and grow spontaneously. This will lead to a multifragmentation of the 
system. We consider here the energy of the system and not the Helmholtz free 
energy because we are describing a post freeze-out situation when we do not have 
a heath bath any more. 

Substituting flT2"|) into flHD and expanding in terms of n\ up to the second 
order, we evaluate the total energy of the perturbed system. For large enough 
systems the terms containing T 2 /fl can also be neglected. Thus the total nuclear 
energy can be written as E = Eq + AE, where 

E E kin (n ) { 4ttVo\ 2 , , 

P + »~ U ° + ^0 ) ( 13 ) 



A A V n 



S 



and A = J d 3 r n(r,r) = J d 3 r no = uqQ. The change of the energy due to the 
formation of the perturbation in the second order of n\ (the first order terms 



cancel due to the mass number conservation: Eq. 12) is then 



AE = AE^- f/3 + -^- J(a+l)(<7 + 2) 7 n£]iV 1 +A£ surf iV 1) 



Ae 



surf 



/i 2 2 

— 5 ==- d r 1 dr 2 ni(ri)ni(r 2 )— — 

N\J | ri — r2 | 



(14) 



where N\ = j d 3 r nf(f). The term Ae sur f is the surface correction due to the 
Yukawa interaction. The sign of AE will determine whether an instability may 
increase or will be damped. 

It is not immediately clear, whether the kinetic energy gives a contribution 
to AE or not. In thermal systems at high temperature and low density, where 
the exact value of the Fermi momentum is not too important we can assume that 
/ = fo + fi, where fi « [n 1 (r, r) /n ] f , and the kinetic energy depends only 
linearly on rii, that is a density perturbation will not cause a change of total 
kinetic energy. For a isotherm, degenerate system, where the kinetic energy is 
nearly proportional to n 5//3 , a perturbation may give a contribution to AE. Here 
we are studying a post freeze-out situation out of thermal equilibrium. Now the 
kinetic energy depends on the form of our perturbed phase space distribution 
function, f s (r,P), we choose or obtain. This may have different characteristics 
depending on the density of our frozen-out system. Thus we will examine both 
situations, the one without the kinetic energy contribution (AEi), and the one 
with (AE 2 ). 

The surface correction, Ae sur f, for the cases of plane wave and spherical 
perturbations considered in Section 2 are: 

4ttV 



Case A 5? 1 - 



Case B (^V /fx 2 ) [1 - g{R*, k)} 

where g is a complicated function of R* and k (see Appendix C). For the used 
forces in Eq. ( |ll|) in the case A AE turns out to be as seen in Ref. |13 . 



^.(, + «.feiMU) Tl<l 

We get a negative change in the energy only for k < fc cr it, where the critical 
value of k is 

^rit - ^ I" AEkm (g+lKg+2) a 

In case B the required negativity of AE can be expressed in a more complicated 
way. 
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HARD EOS, <x=4. 
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Figure 1: The dependence of the energy difference, AE, caused by the formation of 
a bubble of radius R, without the kinetic term on the radius of the bubble for two 
different diffuseness coefficients, a = kTf r /r, and two equations of state. The energy 
difference is evaluated for three post-freeze-out densities no/njy = 0.3, 0.4 and 0.7. The 
maxima of the curves (if any) is at the critical radius, R* 



4. Results 

In the following we present here the results for a spherical droplet perturbation 
considered in Section 2.2. The total energy of the system is given in Eq. (|i~T|) 
and the energy change due to the bubble formation in Eq. (|14|). The surface 
energy, Ae sur f, energy given in Eq. (^) depends strongly on a = krf r /r and on 
the radius i?* of the bubble. To see the effect of the bubble shape on the energy 
change AE we give this change for different a and density values as the function 
of the bubble radius. In Fig. 1 we assumed that the kinetic energy does not give 
any contribution in second order to the energy, in Fig. 2 we considered a ground 
state Fermi kinetic energy contribution. As one can see, the effect of the surface 
term is larger if we have sharper surfaces (larger a values). As one expects, the 
effect of the surface energy, Ae: sur f, decreases for large radii. As a first step we 
considered the solution of Eq. ( |l0l) for n=0, and compared the condition of the 
dynamical instabilities to grow with that of the mechanical instability for infinite 
systems. Without surface term the two condition are the same, as was pointed 



out already in Ref. |p^| . With the surface term the instability region decreases, 
just as in Ref. [13|. In Fig. 3 we give the n - T curve of instability region for 
different a and R* values. With small supercooling first big droplets can nucleate 
and grow, then with stronger supercooling smaller droplets can also be formed. 
For big droplets (R* = 4.5 fm) the effect of the surface term is almost negligible. 
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HARD EOS, 0C=2. 



HARD EOS, cc=4. 
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Figure 2: The same as Fig. 1 with kinetic energy term included 



If we wait longer after the freeze-out in the expanding system, i.e., we increase r 
and thus having a smaller a, (a = KTf r /r) we have the possibility of instabilities 
earlier, with smaller supercooling. The calculations are done both for soft and 
hard equation of state. One sees, that the effect of the surface term is more 
significant for the soft equation of state. 

As the next step we want to determine the break up of the system starting 
from different freeze out densities times and temperatures. A reasonable freeze 
out density should be in the range of no =0.08 - 0.14 fm~ 3 . In our calculation 
we choose nj r = 0.1 fm -3 . Other freeze out densities can be considered simply 
rescaling r/ r , u/ r and Tf r to keep the relation (n /nf r ) 2 ^ 3 = (T/Tf r ) = (rf r /r) 2 
and Hf r /Tf r constant (see the remarks after Eq. (|3])). The freeze out time Tf r 

3 R 2 

defines the flow energy of the system as Eq ow /A = — j- for a system with radius 

5 r fr 

R. We assume, that the total excitation energy of the system is large enough to 



reach the break up densities fl5|| , so if that condition is fulfilled, the freeze out 
time is not defined in the model, and the time scale is not fixed. We followed 
the paths along the trajectories in the n - T plane from different initial freeze 
out configurations. The instability condition Eq. ( |TDj ) is examined along the 
trajectories, and the time of the solution for growing instability (break up) can 

Ifr 

n 

In Fig. 4 we show the part of the trajectories of the post-freeze-out expansion 
denoted by dotted line where the instability condition is fulfilled and the energy 
change, AE ly is negative. We found, that whenever condition ( ]T0| ) is fulfilled, this 



be expressed as At = r — Tf r = Tf r 



- 1 



after the freeze out. 
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Instability region 
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Figure 3: The dynamical and mechanical instability regions. The region of mechanical 
instability is bordered by the isothermal spinodal, where the isotherm sound speed 
vanishes and becomes imaginary. For large size (radius R) of the critical bubbles or 
for systems without surface energy the dynamical instability region is the same as the 
mechanical one. The break up instability depends on the radius R and on the stretching 
of the system parameterized by a 



energy change is always negative. The more strict condition, using AE 2 , however, 
excludes some part of the trajectories (solid line) at temperatures above 10 MeV 
for the hard equation of state. For the soft equation of state there is no such 
exclusion, nevertheless, the instability region is smaller. We give the results for 
the two equations of state and different parameters for both. For comparison we 
show the boundary of the isothermal mechanical instability (which corresponds 
to the R* — > oo situation). In the parameter region we count as physical (R* ~ 
2-3 fm, k pa 4 - 6 fm -1 ) there are no significant changes on this parameters. 

Following the trajectory of the post-freeze-out expansion starting from a given 
initial nf r ,Tf r configuration the instability condition Eq. fllCf ) has solutions for 
different S or k values. For high densities and temperatures S is negative, that 
is there are no real solution for k. As the density (and correspondingly the tem- 
perature) decreases it continuously becomes positive. The speed of growth of the 

instabilities is determined by k = -JlT^l J — / r . Evaluating this expression we 

r r y m{R*Y 

assumed that r^/r ph 1. In Fig. 5 k is shown along given expansion trajectories 
(fc=4 fm" 1 R= 2 fm, T fr = 1 MeV, 16 MeV and 21 MeV) for the hard equation 
of state. The speed k, the instabilities are developing with, is changing, first it 
increases as the nucleus evolves to smaller densities, later it decreases back to 
zero. The time scale of the expansion for the radius of the perturbations from 
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Figure 4: The trajectories of an adiabatic post-freeze-out expansion starting at the 
same, rif r = 0.1 fm™ 3 freeze out density. The different lines are originated from different 
freeze out temperatures. The solid lines correspond to the case where the energy 
difference with the kinetic term is negative, the dotted line is the case where the 
energy difference without the kinetic term is negative. The latter ones define the wider 
region. The possibility of dynamical post-freeze-out instability opens only after some 
penetration into the domain of the isothermal spinodal, while thermally dominated 
homogeneous nucleation may start immediately, although slowly 
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Figure 5: The growth parameter, k, along trajectories of post-freeze-out expansion 
starting at nj r = 0.1 fm -3 and at the given freeze-out temperatures 



Fig. 5 gives m 20 fm/c for the k ~ 0.06. However, the growth rate of the insta- 
bility from Eq. (|j) is a double exponential: n c exp (kTfr/rth R*e KAiT ~ Tth) ), which 
is much faster. This region of the exponentially developing instabilities breaks 
up the system. 



Appendix A: Dispersion relation for plane wave perturba- 
tion 



Let us consider a plane wave perturbation in another form than in Ref. [JT2 
emphasizing our time scale 



/i(r,r,p) = f s (r,P) exp 



% — —r + k(t - Tth) 

T 



(15) 



where it is taken into account, that the wave number, k, scales with r, as all other 
parameters of the flow. Thus A; is a constant, independent of r. Here fi is a phase 
space distribution which should satisfy the Vlasov equation @, and it represents 
a plane wave with growing amplitude if k is a real positive number. Since we 
assume small perturbations only, this solution is valid only for a short time after 
the instability starts to grow. (Frequently in similar studies the perturbation is 
studied in a form containing an exp(ia;r) term, and the opening of the channel of 
instability is indicated when u becomes imaginary. Thus perturbations preceding 
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r t h will not grow. This approach is basically equivalent to the one presented here, 
however, we wanted to emphasize that the instability may grow only after r^.) 

As it was already mentioned we also assume that / s (r, P) has a characteristic 
time dependence on the slow scale (i) of the post freeze-out expansion, so that 
its time-derivatives are negligible compared to other time-derivatives of /i(r, f,p) 
corresponding to rapid inequilibrium processes (ii) like exp[/?(r — r t h)]- Further- 
more, we assume that f s depends on the LR momentum, P only, i.e., it does not 
have any other dependence on f other than what is included in P. 

We search for such solutions of the Vlasov equation, fi(r,r,p) and f s (r,P). 
We assume that such solutions can be obtained only some time, r, after the 
freeze-out at rj r , i.e. at r th > r/y. 

Using the above form, (|5|), of perturbation the density change arising from 
this perturbation is: 



ni[r,r 



n s (r) exp[i — — r + k(t 



Tth)]- 



Inserting a plane wave perturbation Eq. (|5|) into the Vlasov equation we obtain 
for the force used in Eq. fllTD 



f s (r,P) 



ikTfj. . _ 

K H — [p 



mr 



mf/r) 



dU(k) AkTf r ( / n 2 \ [p-mr/r] 



[p-mf/r] 2 Heff 



2mT f 



eff 



where 



dU(k) 



-2/3 + 7 ((7 + l)((7 + 2) n% 



T eff 

AttVo 



(16) 



as in Ref. |TJ. 

Note that the momentum appears only inside expressions of the LR momen- 
tum, p — mf/r, thus we can integrate it out in the LR frame. Thus from (|16|) 
we can express f s , and integrating it over the LR momentum, P = p — mf/r, we 
obtain n s (r), which then can be eliminated from both sides yielding: 



1 



dU(k) 

T eff C dn 



d 3 PfZ(P)- 



%Tf r kP jr 



inn 



ikrf r 



P 



exp 



P 2 



2mT, 



eff 



eff 



(17) 



We can separate the variable of the integral to a parallel, P\\, and an orthogo- 
nal, P±, component with respect to k, and the integration over the components 
perpendicular to k can be performed. This yields 



1 



-2TtmC 



dU(k) 
dn 



dPi 



ikPuTf r /T 
+ ikP\\Tf r /r 



1 + exp 



Veff 



2mT, 



eff 



T, 



eff] 



18) 
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Only symmetric functions contribute to this integral, so we symmetrize it by 
multiplying both the numerator and denominator by ttik — ikP\\Tf r /r, and then 
dropping the antisymmetric term we end up having 



dU ' 



f r \\ 



^11 _ MS 
2mT eff T eff 



• (19) 



Introducing a new variable, y = Pu/(2mT e ff), a straightforward calculation will 
lead to the dispersion relation in the integral form: 

OO 

1 = -2ixmCj2mT eff ^- [ - ^ V (20) 

v ;/ dn J (y + s)[l + ey-^/Tfry v 1 

where s = mK 2 r 2 /(2k 2 r] r T eff ) = mK 2 r 4 /{2k 2 T fr Tf r ) . 



Appendix B: Dispersion relation for spherical perturbation 



For the sake of simplicity let us first study a spherical cusp perturbation centered 
around some interior point r c (r) = r r/rj r of the type 



fi = exp [-kr fr \f- r T/T fr \/r + k(t - r th )] f s , 



(21) 



where the center of the perturbation moves along the scaling expansion, and this 
center was at point r at the time of the freeze-out. As we discussed it in the 
introduction we can assume that r = 0, so without loosing the generality of the 
assumption, \ f— r | — > |rj = r, since the interior points of the expanding nuclear 
system are equivalent. Thus 



fi = exp [-kr fr r/r + k(t - T th )\ f s , 



(22) 



can serve to study the perturbation just as well. Although this functional form 
of perturbation has a singularity in the center this will not be essential for our 
study, and it could be removed by assuming more complicated functional forms 
for the perturbation which would not show such a singularity. 

However, including the center of the perturbation, r c , explicitly will allow 
us later to discuss interactions (e.g. fusion, repulsion, etc.) of two (or more) 
elementary perturbations. Thus we will follow this somewhat more complicated 
derivation although it is not necessary at the moment. 

We will not explicitly define the form of f s at this stage, unlike in the case 
of plane wave perturbations. Instead we will consider the integrals of fi and / s . 
First the norms: 



= J d 3 p h 



exp 



T 



T 



r — r - 



+ K(T ~ T th ) 
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and n s = J d 3 pf s , where we require that n s = n s (r) does not depend on r. Second 
the projections orthogonal to (r — f c ), i.e., 



gi{r,r 



;P||) = 2it j P±dP ± f 1 =exp 



T 



r - 



+ KyT - Tth) 



Here g s (r,Pu) = 2n J P±dP±f s , where we require that g s = g s (r,Pu) does not 
depend on r. We chose our coordinate system so that Pi is parallel to (r — r c ), 
and Pj_ is orthogonal to it. These constraints are the required implicit constraints 
on the choice of f s . 

Inserting Eq. (^l|) into the Vlasov equation and integrating it over d 2 P± we 
obtain 



kTf, 



Pi 



f r/ Tf r 



kr 



vtlt r 



f r/ Tfr 



i^T /Tf r 



2*cn s d 4- ^> ,:~ r } r ( r/ - rPu 

on t \r — r T/Tf r \ 



1 + exp 



kfp f-fpT/Tfr 
T \r-foT/Tf r \_ 

-ix — 1 



2mT, 



eff 



T eff\ 



0. 



(23) 



Performing products in the first term and taking into account that P ,t r ° T / T f r 

■ \r— rQT/Tf r \ 



where for an averaged Yukawa surface term U is only the function of n. 

Performing ] 
= Pi I we obtain 

kr frr ^ , f — r r/T f 



9s 



'fr 

rriT 



[p — mf/ t] ■ 



f r/ Tfr 



K + 



kTfr n \r- r T/r fr 
1 -P\ 



so that Eq. 



mr " \r — for/r^v 
takes the form 



9a 



K+— P \\ 
VTLT 



„ dU kTf r { 

27rCn s - } - Pnh + exp 

on t 



K + 



Pi 



kT 



fr 



TTIT 



P 



2mT ( 



eff 



T, 



eff 



-1 



0. (24) 



We see that the position of the center of the perturbation, f , has dropped 
out. We indicated this symmetry already in the introduction when we pointed 
out that in the spherical scaling expansion all interior points are equivalent in 
the sense of their LR features. We would obviously get the same result assuming 
f = from Eq. ( Ell ) on in the course of this derivation. 

Dividing both sides by (k + -^Py) and integrating over dP\\ leads to 



2^ *?> 

on t 



oo 


kTf r 


1 \ 


jP\\dP\\ 


< 1 + exp 


-co 


rriT 



flfr 



-ix — 1 



2mT, 



eff 



T 



fr] 



(25) 
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Multiplying both the denominator and the numerator by k — ~^-P\\ and then 
dropping the antisymmetric part, which does not contribute to the integral we 
obtain 



au 7 
-AnmC— dP\\Pfi 
onJ 11 
o 



m 2 T 2 n 2 
k 2 r 2 r 



-Ml 



1 + exp 



II 



2mT f 



eff 



T 



(26) 



The same way as we got the dispersion relation in the case of plane wave pertur- 
bation from Eq. fll9"f ) we get now the relation 



dU 



2irmC\/2mT e ff — — 



dy y/y 



On J (y- s) [1 + e y ~Vfr/T Sr ^ ' 



(27) 



where s = mK 2 r 2 / \2k 2 T 2 r T eff ) = mK 2 r 4 / (2k 2 T fr rj r ) . 

If we have a spherical droplet perturbation of a finite radius described by 
Eq. (H) instead of a cusp, Eq. fl2~ID , the dispersion relation can be obtained from 
( |27| ) by making the transformation k — > K,kR*Tf r /T t h arising from comparing 
the form of the two perturbations ([H|) and (||). This leads to exactly the same 
equation as the equation above (|2"7|) for the spherical cusp perturbations, except 



that in place of s we have 5* = mK 2 (R*) 2 T 4: /(2Tf r Tj r T 2 h ) in the expression. 



Appendix C: Spherical droplet with Yukawa forces 

For the profile (|9|) the Yukawa force can be written as follows 



V Y uk(r) 



4ttV 



-e~» R sinh(//r) f — - 



R 



pL + a /i 2 (fx + a 



for r < R(t), where a = krf r /r, and 



2a 



a 2 — jj 2 r(a 2 — fx 2 ) 2 



e («- M )(r-H) (R I R 1 
+ ~ ~ + + 



+ 



2r/i 



[i [i* a — ii 
'R 1 R 



a — id) 2 j 
1 



2rfie( a+ ri R \fi fi 2 a + /i (a + /j) 2 



for r > R(t). 

The Yukawa energy can be written as 



3 



Yuk 



d^r V Yv k{r)ni(r) 



AWq 2 47T 



2 "°3 



3 3 2 a 2 -a/i-/i 2 

it ri — 

2 jia{jJL + a) 



(28) 



(29) 
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+ 



1 R _2_L_ 

2 \ fi{fi + a) 

3 / (2a + fj,)fj, 
2 V2a 3 (/i + a) 2 



-2/xiJ 



Ra 




(30) 



j(i(/i + a) (/i + a 



One sees that substituting expression fl28|) and ( p0|) into the Vlasov equation, 
the perturbation (|T0|) is not a solution of it. However, for sharply decreasing 
surfaces aR > 1 and we can consider instead of Vynk the average of it. That is, 
we use a surface term, which is given as the average of the Yukawa interaction. 
Instead of ([3(]) we introduce for r > R{r) (see Eq. (|TTD) 

4ttVo 



11 



surf 



AttV 



a 



1 — V^Yuk 
R + 



(31) 



nnr) 



n{n + a) R 2 + ^ + ^ 



The surface energy term in Eq. ([11]) can be written using fl3~0|) and neglecting 
the terms proportional to e~ 2fJjR (these terms are small) as 



AirVr 







E 



Yuk 



surf 



1 4vrV 



R 2 a 



+ 



R 2a + fi 4a + fi 



a 



fi 3 R jj, + a 



2/j(n + a) ' (/i + a) 2 
for aR > 1 . 



+ 



2/i 4a(/i + a) 2 2/i 2 



(32) 



Appendix D: The mechanical instability region 

The total energy density of the system can be written as e pot {n) + ep, with the 
potential energy density e pot (n) and the kinetic (Fermi) energy density cf- The 
latter can be written as cf = const. T 5//2 F 3 / 2 (/i/T), using the integrals Fi/ 2 (n) = 

/ dx 1+e *p( x _ v ) ■ The density can be expressed as n = const. T 3 ^ 2 Fi/ 2 (fx/T). 



Isotherm expansion, dT = 0, leads to the change of Fermi energy 

dep 



3T^pm dn 



The pressure should be calculated from the free energy density f(T,n) = e — Ts, 
with the entropy density s = \cf/T — nfi/T: 

2 df/n dc nnt {n) 



p = n 



dn 



n- 



^pot 

dn 



ipot{n) + -e F 
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The region of mechanical instability where the derivative of the pressure above 
with respect to the density at constant temperature is negative: 

dp = d 2 e pot (n) F x/2 {n/T) _ dU F 1/2 Qm/T) 
dn dn> + F. 1/2 (jjl/T) dn + F_ 1/2 (/i/T) 

The onset of the instability is determined then by 

1 = -2 cT ' a^ F -^m , ° = &{w) (33) 

The condition of the dynamical instability derived in Appendix B is 

= 27rr /r m g [^fr^L 7 ^v 7 ^ 

r (2tt^)3 V f r dnJ (y- s) (1 + exp(y - fi eff /T eff )) 

1 /2 

substituting Tf r /r = (T/Tyy) 1 for the critical mode (s = growing) one gets 

1 /2m\ 3 / 2 g ^/odU „ , , m . 

1 = -2(if) ^ftr'-v»0w//^/) • (34) 

which is the same as the condition for the isothermal spinodal. 



Appendix E: Critical radius 

Let us introduce the notation (nf) 4^ R 3 = N 1 and rewrite Eq. (02] ) as 

Aiw = ^-i^-^ = .MJf. (35) 
\i A Kfi{n + a) 

Using Eq. ( [HQ the total energy change due to the formation of a droplet of size 
R can be cast in the form 

AE(R) = A£ lln " (f + ^ " (g + 1 f + 2) -K) M"?> * + ^ 

=6 

Thus, assuming that the contribution of kinetic energy is independent of R we 
obtain the critical radius from dAE/dR = 



At _R* rit the function AE(R) has its maximum, thus bubbles or droplets smaller 
than R* vit shrink while larger than R* rit grow. 
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